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ABSTRACT 

In this paper we report on PD-SPH the new tree-sph code developed in Padua. The 
main features of the code are described and the results of a new and independent series 
of 1-D and 3-D tests are shown. The paper is mainly dedicated to the presentation of 
the code and to the critical discussion of its performances. In particular great attention 
is devoted to the convergency analysis. The code is highly adaptive in space and time 
by means of individual smoothing lengths and individual time steps. At present it 
contains both dark and baryonic matter, this latter in form of gas and stars, cooling, 
thermal conduction, star formation, and feed-back from Type I and II supernovae, 
stellar winds, and ultraviolet flux from massive stars, and finally chemical enrichment. 
New cooling rates that depend on the metal abundance of the interstellar medium are 
employed, and the differences with respect to the standard ones are outlined. Finally, 
we show the simulation of the dynamical and chemical evolution of a disk-like galaxy 
with and without feed-back. The code is suitably designed to study in a global fashion 
the problem of formation and evolution of elliptical galaxies, and in particular to feed 
a spectro-photometric code from which the integrated spectra, magnitudes, and colors 
(together with their spatial gradients) can be derived. 
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The conventional picture of star formation in elliptical galax- 
ies is one in which the galaxies and their stellar content 
formed early in the universe and have evolved quiescently 
ever since. This view is supported by the apparent unifor- 
mity of elliptical galaxies in photometric and chemical ap- 
pearance (cf. Matteucci 1997 for a recent review) and the 
existence of scaling relations, e.g. the fundamental plane (cf. 
Bender 1997). In contrast, the close scrutiny of nearby ellip- 
tical galaxies makes evident a large variety of morphological 
and kinematic peculiarities and occurrence of star formation 
in a recent past (Schweizer et al. 1990, Schweizer & Seitzer 
1992, Rampazzo et al. 1996). All this leads to a different 
picture in which elliptical galaxies are formed by mergers 
and/or accretion of smaller units over a time scale compa- 
rable to the Hubble time. Furthermore, strong evolution in 
the population of early type galaxies has been reported by 
Kauffmann et al. (1996) which has been considered to sup- 
port the hierarchical galaxy formation model (Kauffmann et 
al. 1993, Baugh et al. 1996). 

Tracing back the formation mechanism of elliptical 
galaxies from the bulk of their chemo-spectro-photometric 
properties is a cumbersome affair because studies of stel- 



lar populations in integrated light reveal only luminosity 
weighted ages and metallicities, and are ultimately unable to 
distinguish between episodic (perhaps recurrent) and mono- 
lithic histories of star formation and between star formation 
histories in isolation or interaction. 

As nowadays most of the properties of elliptical galax- 
ies have been studied with sophisticated chemo-spectro- 
photometric models in which the dynamical process of 
galaxy formation is reduced to assuming either the closed 
box or infall scheme and a suitable law of star formation 
(Arimoto & Yoshii 1987, 1989; Bruzual & Chariot 1993; 
Bressan et al. 1994; Worthey 1994; Einsel et al. 1995; Tan- 
talo et al. 1996; Bressan et al. 1996; Gibson 1996a,b,c; Gib- 
son 1997; Gibson & Matteucci 1996; Tantalo et al. 1997). 

In contrast, the highly sophisticated dynamical mod- 
els of galaxy formation (Hernquist & Katz 1989; Davis et 
al. 1992; Katz 1992; Steinmetz & Miiller 1993; Navarro & 
White 1993; Nelson & Papaloizou 1994; Katz, Weinberg & 
Hernquist 1995; Navarro & Steinmetz 1997; Haehnelt et al. 
1996a,b; Steinmetz & Mueller 1995; Navarro et al. 1996; 
Steinmetz 1996a,b and references), owing to the complexity 
of the problem are still somewhat unable to make detailed 
predictions about the chemo-spectro-photometric properties 
of the galaxies in question. 
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Attempts to bridge the two aspects of the same prob- 
lem are by Theis et al. (1992 and references) and most re- 
cently by Contardo, Steinmetz & Fritze-v-Alvensleben (1997 
in preparation). 

This paper is the first of a series dedicated to the study 
of this complex problem in a self-consistent fashion, in which 
the formation and evolution of galaxies, together with their 
chemo-spectro-photometric properties, stem from a unique 
model able to predict a number of observable quantities to 
be compared with real data. The project aims to build dy- 
namical models of disk-like and elliptical galaxies by means 
of the SPH (Smoothed Particle Hydrodynamics) technique, 
in which the effects of different initial conditions and ba- 
sic physical processes, such as star formation, heating of 
the gas by various mechanisms (supernova explosions, stel- 
lar winds, UV fluxes), cooling of this by radiative atomic 
and molecular agents, interplay between dynamics and ther- 
modynamics (feed-back), and chemical enrichment of the 
inter-stellar medium are taken into account. The models al- 
low for the presence of dark matter and the effect of this 
on the gravitational field, which is described using Barnes 
& Hut (1986) treecode. The rates of star formation and 
chemical enrichment as function of space and time result- 
ing from the dynamical models are meant to feed chemical- 
spectro-photometric models which should provide the kind 
of data (luminosities, integrated spectral energy distribu- 
tions, broad band colors, line strength indices, etc.) to be 
compared with real observations for nearby and high red- 
shift galaxies. 

In this paper we present the first step toward this ar- 
ticulated analysis of the problem, i.e. the tool to construct 
dynamical models, and describe in some detail the numerical 
algorithm and the key physical ingredients. The dynamical 
models are based on the SPH technique, one of the most 
popular and efficient tool of modern studies of numerical 
hydrodynamics, and in astrophysical context of structures 
and galaxy formation. In the following we will refer to the 
SPH code we have developed as the Padua SPH (PD-SPH). 

When SPH is coupled with an efficient scheme to com- 
pute gravity (like tree codes), it becomes a powerful tool to 
investigate physical situations of enormous complexity such 
as for instance, the formation and evolution of a galaxy (in 
the isolation picture) and/or the interaction- merge of galax- 
ies (in the hierarchical scenario). The lagrangian nature of 
the method allows us to follow the evolution of physical 
quantities at all scales, provided that sufficient resolution 
is ensured. 

The plan of the paper is as follows. Section 2 gives a 
detailed description of the code. Section 3 discusses tests in 
one dimension, i.e. the results for the Riemann problem and 
the convergency analysis, whereas section 4 contains tests 
performed in three dimensions. Section 5 presents the first 
astrophysical application, i.e. the collapse of a galaxy made 
of gas and dark matter in presence of initial solid body rota- 
tion. Section 6 thoroughly presents non-adiabatic processes, 
such as thermal conduction, radiative cooling, and heating 
by Type I and II supernova explosions, and stellar winds 
and ultraviolet radiation from massive stars that have used 
as source or sink of energy in our models. We adopt cooling 
rates that take the effect of metallicity into account. In Sec- 
tion 7 we describe the results obtained from our cooling pre- 
scription and compare them with those from metal indepen- 



dent cooling rates. Section 8 describes our prescriptions for 
the star formation rate and feed-back (rate of energy input 
from supernovae, stellar wind, and UV emission), and rate of 
chemical enrichment. Section 9 presents the results for two 
disk-like galaxies: the first case is with no feed-back of any 
type, whereas the second case is with feed-back. These mod- 
els are not meant to represent real galaxies, first because the 
initial conditions are not derived from a cosmological con- 
text, second because the total number of particles is small 
owing to limitations in the computing facilities. The models 
are indeed meant to test the code response to various physi- 
cal inputs. Nevertheless, despite the present limitations, the 
results are very promising. Finally, Section 10 presents some 
concluding remarks and outline the perspectives of future 
studies. 



2 THE PD-SPH CODE 

PD-SPH is a TREE-SPH code, written in Fortran 90 (cf. 
Carraro 1995 & Lia 1996), conceptually similar to those 
described by Hernquist & Katz (1989) and Nelson & Pa- 
paloizou (1994). A 1-D version of the SPH code is in par- 
allel form, and it has already been tested on the Cray T3D 
parallel computer hosted by CINECA (Beninca & Carraro 
1995). 

Schematically, PD-SPH uses SPH to solve the equations 
of motion for the gas component, and the Barnes & Hut 
(1986) octo-tree to compute the gravitational interactions. 
In this section we describe in some details the structure of 
the code and its basic ingredients. 

It is widely known that SPH is a method standing on 
two basic ingredients: an interpolation using a kernel, and 
a Monte-Carlo evaluation of the physical quantities (Mon- 
aghan 1985). All relevant variables are interpolated in the 
following way: 

</(f)> = Jw(r-?,h)f(/)dS, (1) 

where the integral is extended over the whole space and 
IF is a function generally referred to as the interpolating 
kernel; h is the smoothing length determining the spatial 
region within which variables are smoothed, and governing 
the spatial resolution of the method. 

The kernel is a spherically symmetric function strongly 
peaked at r — f . It is easy to prove (Benz 1990) that these 
properties make SPH a second order technique, in which the 
estimate of any function is given by 

(m)=f(7)+c(V 2 f)h 2 + 0(h 3 ). (2) 

This estimate is made in lagrangian formalism summing up 
the contributions of all the particles within 2 x h, which 
represent the so-called neighbors. 



2.1 Interpolating kernel 

Three types of kernel can basically be found in literature: 
exponential, gaussian and spline; see Benz (1990) for more 
details. We have adopted the spline-kernel of Monaghan & 
Lattanzio, (1985): 
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W(r,h) 




3 u> + \u\ 



se < u < 1; 
se 1 < u < 2; 
otherwise, 



(3) 



and the supergaussian-kernel of Gingold 
(1982): 



W{r,h) 



1 



7T 3 / 2 /l 3V 2 



(- — u)exp(—u ), 



Monaghan, 



(4) 



where « = | //*. 

The super-gaussian kernel is built up forcing the sec- 
ond order term in eq. (2) to be zero, It can be shown that 
this kernel corresponds to a fourth ord er in terpolation. This 
kernel becomes negative for r/h > y/5/2 (Benz 1990) and 
does not find an immediate physical interpretation. 

Moreover, while the spline kernel is defined on a com- 
pact support, the super-gaussian kernel needs an artificial 
cut-off to limit the number of neighbors utilized in the es- 
timate of the physical quantities. First, we checked both 
kernels (see Section 3 for details) and finally adopted the 
spline-kernel. In our code it is stored as a look-up table. 



2.2 Evolution of the smoothing length and the 
search of neighbor particles 

In SPH the spatial resolution is fixed by the smoothing 
length h. In our code h is kept variable both in space and 
time according to Benz (1990) 

dh I K dpi 

^ = ~3^ = 3 WW - (5) 

This equation is added to the set of equations to be solved at 
any time-step. The main drawback of this approach is that 
in situations in which strong density gradients are present, 
it is not possible to keep under control the number of neigh- 
bors. Typically this number should be about 40-50 (Stein- 
metz & Muller 1993). 

One possible solution is to artificially increase or de- 
crease h in order to keep fixed the neighbor number, but 
this involves several iterations if - as in our case - the tree- 
code is used to look for neighbors. As a consequence of this, 
the computational time gets often unacceptably long. 

To cope with this difficulty, we have decided to maintain 
the above scheme when the neighbor number is lower than 
40, and switch to the Nelson & Papaloizou (1994) formalism 
when the neighbor number is greater than 80. Accordingly 
the new smoothing length h is 



hi = 



1 



N far ^ 2 



(V f ar 

^ 1,- 

- \ri — r„ 



(6) 



where Nf ar are the n most distant nearest neighbors. We 
assume Nf ar — 10. 

There is another drawback in the Benz (1990) formal- 
ism that deserves some attention. The space-time variation 
of the smoothing length generally hampers the strict conser- 
vation of energy when terms involving the space derivatives 
of h are not included. V/i terms are essentially small cor- 
rection (cf. the thorough discussion by Nelson & Papaloizou 
1994), which do not significantly improve upon the conser- 
vation of energy. 

Moreover the momentum conservation has been secured 
at an acceptable level of confidence also by symmetrizing the 



search of the nearest neighbors. In brief, a candidate particle 
i is considered as a neighbor of j if the following conditions 
are met 

\n - fj\ < 2 • hi 
or 

\n - f*,- < 2 • hj. 

2.3 Hydrodynamical equations 

In PD-SPH, the particle forces and the specific energy are 
computed by means of the following equations 



dvi 
~dt 

1 

2 

and 



JV 



PiPi 



ViW{n 3 ,K) +V i W(r ih h j ))) 

Vp^ i tt \. 



(7) 




Pi Pi 



\ (ViW( rij , hi) + ViWir-rj.hj)) + r R - 



(8) 



Ac 

1 p 

in which we adopt the arithmetic mean for the pressure gra- 
dient as in Hernquist & Katz (1989). Furthermore, in eq. 
(^) the first term represents the heating rate of mechani- 
cal nature (it is shortly indicated as Fm), the second term 
Fr is the total heating rate from all sources apart from 
the mechanical ones, and the third term Ac/p is the to- 
tal cooling rate by many physical agents. These latter two 
will be discussed below in more detail. In the above equa- 
tions Vij = Vi — Vj, and Ylij is the viscosity tensor defined 



-aCijUij+Pti,- 



where 



if (Hi — Vj) ■ (n — fj) > 
otherwise, 



p*ij 



hij(vi - Vj)(ri - r f ) 



(9) 



(10) 



Here dj = 0.5(c;+Cj) is the sound speed, htj = 0.5(hi +hj), 
and a and (3 are the viscosity parameters, usually set to 1.0 
and 2.0, respectively. The factor e is fixed to 0.01 and is 
meant to avoid divergencies. 

As amply discussed by Navarro & Steinmetz (1997) this 
formulation has the disadvantage of not vanishing in the 
case of shear dominated flows, when V ■ v = and V x 
v ^ 0. In such a case, spurious shear viscosity can develop, 
mainly in simulations involving a small number of particles. 
To reduce the shear component we adopt the Balsara (1995) 
formulation of the viscosity tensor 

fi + Si 



n = Uu x 

z 

where /; is a suitable function defined as 
< V ■ v >i I 



(11) 



(12) 



Si = — = — 

<V-i7>i| + |<VxiT>i|+ rjd/hi 

and rj ~ 10~ 4 is a parameter meant to prevent numerical 
divergencies. 
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Figure 1. Result of the shock tube problem using super-gaussian kernel (left panel) and spline kernel (right panel). In both cases 400 
particles have been used. The arrow indicates the so-called blip in the pressure profile 



2.4 Time integration and stability criteria 

Particle positions and velocities are updated, as in Hern- 
quist & Katz (1989), using the leapfrog integrator and the 
multiple time-step scheme. The integration is of the second 
order in time, and proceeds as in the classical scheme. 

Firstly an estimate of the velocity is obtained 



from 



-ra+l/2 
Vj 



= t£* +0.5At;a™- 1/2 - 



(13) 



This is used to compute time-centered accelerations, a" +1 ^ 2 , 
from which particle velocities and positions are updated 



v" +1 =i% + AUa" +1/2 



^+1/2 = 



(14) 
(15) 



The energy equation is explicitly solved, unless sink or 
source terms are present, and particle energies are advanced 
in the same manner as positions. 

Time steps are calculated according to the Courant con- 
dition 
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Figure 2. The number of particles involved in the blip versus the 
total number of particles in the simulation. Nt,up is on logarithmic 
scale 



: C- 



hi 



hi V • w + Ci + 1.2(aa + (3maxj \fj,ij\) ' 



(16) 



with C w 0.3. 

In presence of gravity, a more stringent condition on 
the time steps is required. According to Katz, Weinberg & 
Hernquist (1995), the additional criterion has to be satisfied 



\v\ a 



(17) 



where e is the gravitational softening parameter and r\ is 
another parameter usually set to 0.5. The final time step to 
be adopted is the smallest of the two 



AU = MIN(Atc,i,AtG,i) 



(18) 



While the use of the multiple time step scheme may signif- 
icantly speed up the code (Steinmetz 1995), sometimes the 
interpolations of non active particles may affect the energy 
conservation (see Section 4). 



3 1-D NUMERICAL TESTS 

The classical 1-D test to which SPH codes are compared 
is the Riemann shock tube problem. We use this test not 
only to check whether our code works properly, but also to 
analyse its convergency performances. The classical refer- 
ence for this test is Gingold & Monaghan (1983), to whom 
we refer for the initial conditions of the problem. The time- 
steps and the smoothing lengths are constant, and fixed to 
0.05 and 0.025, respectively. In other words, the smoothing 
lengths are set to be two times the maximum inter-particles 
separation. The particle masses (0.003125 in our case) are 
obviously chosen in such a way that the density and pressure 
profiles are matched (Gingold & Monaghan, 1983). The ker- 
nel is suitably normalized using the normalization constant 
(see for instance Hernquist & Katz 1989.) 



Figure 3. Blow-up of the contact discontinuity in the pressure 
profile. The left panel (a) shows the results obtained from adopt- 
ing the super-gaussian kernel and the force term from Gingold & 
Monaghan (1983). The central (b) and right (c) panels show the 
same but for standard formalism for the pressure gradient and the 
super-gaussian and spline kernel, respectively. The dashed line is 
the analytical solution 
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Figure 4. Convergency test. The differences A between the nu- 
merical and analytical solutions for pressure and velocity profiles 
as a function of the total number of particle in the simulation 

In Fig. |j] we show the results for pressure, density, veloc- 
ity, and entropy of a simulation of the shock tube problem 
obtained using 400 particles and both the spline and super- 
gaussian kernel (right and left panels, respectively). The aim 
is to check whether using different kernels it might possible 
to eliminate the so called blip in the pressure profile (indi- 
cated by an arrow in Fig. Q). The nature of the blip is well 
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text. Open circles and crosses are for the minimum and maximum pressure. The continuous line represents the analytical solution 



known, and it is the signature of the contact discontinuity 
in the shock. In practice the pressure is discontinuous at 
this point, and the correspondent smoothed value < Pr > 
is expected to be somewhat lower and higher at the left and 
right side of the contact discontinuity, respectively (Gingold 
& Monaghan 1983). Even if in our simulation the shock is 
broadened over a 4 x h region, the blip occurs as the signature 
of the unphysical step in the pressure. 

We find that even using the super-gaussian kernel the 
blip albeit smoothed out does not disappear. In contrast to 
Gingold and Monaghan (1983), we are not able to eliminate 
the blip even using their new formulation for the pressure 
force 
VPr 



p 7-2 [V(Prp) + (7 - l)PrVp] . 



(19) 



Instead of trying other expressions for the pressure gradient, 
once adopted a certain combination of kernel and pressure 
gradient, we prefer to run simulations with a different num- 
ber of particles to check whether the SPH code can eventu- 
ally recover the analytical solutions when larger and larger 
numbers of particles are used. 



In Fig. ^ we show the number, Nbu p , of particles in- 
volved in the blip as a function of the total number of parti- 
cles N to t in the simulation. Nbu v is the number of particles 
counted in a searching sphere of radius equal to 2 x h near 
the blip location, and it is normalized to the total number 
Ntot. The results are as expected, in the sense that Nuip 
strongly decreases at increasing Ntot, or equivalently at de- 
creasing smoothing scale. These experiments show that the 
blip survives because of the still insufficient resolution. 

In Fig. ^ we show the results of simulations performed 
using different combinations of kernels and gradient pres- 
sure terms. The region of interest is blown up for the sake of 
better understanding. Panel (a) is for super-gaussian kernel 
and the functional expression for by Gingold & Mon- 
aghan (1983). Panels (b) and (c) show the same but using 
the standard formalism for the pressure gradient, and adopt- 
ing super-gaussian and spline kernel, respectively. In all the 
three panels the blip is clearly visible. 

We notice, in particular, that the first combination of 
kernel and pressure gradient yields the best results, even if 
the numerical solution strongly departs from the analytical 
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Figure 6. Adiabatic collapse: the evolution of the radial velocity profile (in units of u* = (^r) 1 / 2 ) as a function of the age. The selected 
ages are the same as in Hernquist <fc Katz (1989) 



one. In the case with super-gaussian kernel and standard 
pressure gradient (panel b), the step in pressure is exactly 
the opposite of what is expected. Finally, in contrast with 
Gingold & Monaghan (1983), the inclusion of thermal con- 
duction (cf. Lia 1996) does not eliminate the pressure step. 
Most likely, the blip is a feature intrinsic to the SPH smooth- 
ing itself. 

As far as testing the convergency is concerned, we pro- 
ceed as follows. We select four regions in the pressure profile, 
derive the maximum and minimum values from the numeri- 
cal solutions at increasing number of particles in the simula- 



tion, and compare them with the analytical one. We choose 
the pressure as test physical quantity because of its sensi- 
tivity to the numerical technique. The pressure in fact is 
computed from the smoothed density and energy, instead of 
being directly smoothed out. 

The four selected regions are the undisturbed warm 
fluid at the left of the shock (1), two regions near the shock, 
to the left (2) and to the right (3), respectively, and a region 
just to the right of the blip (4). These regions are centered 
at the linear coordinate X — -0.6, -0.2, +0.04 and +0.2, 
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Figure 7. Adiabatic collapse: the evolution of the specific internal energy profile (in units of m* = — ^-) as a function of the age. The 
selected ages are the same as in Hernquist & Katz (1989) 



respectively. The particles considered in the pressure evalu- 
ation are those contained in the 2 x h searching sphere. 

The differences between numerical and analytical solu- 
tions are shown in Figs. 4 and 5. 

In Fig. ^|we show the pressure (left) and velocity (right) 
differences A between numerical and analytical solutions at 
increasing number of particles in the simulation. The plot- 
ted differences are the mean of their values in the four test 
regions. As expected, the largest differences occur in the 
pressure profile. The four panels of Fig. ^ show the maxi- 
mum and minimum values of the pressure in the four regions 



at increasing total number of particles as indicated. In each 
panel the horizontal line shows the analytical case by defi- 
nition. 

All these numerical tests prove the quick convergency of 
the numerical solution to the analytical one, and compared 
with the similar experiments by other authors, show that 
satisfactory agreement is achieved. 

Finally, we like to report that the strong double-shocks 
tube experiment (cf. Steinmetz & Miiller 1993) has also been 
successfully performed (Carraro 1995). 
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3-D NUMERICAL TESTS 



4.1 Gravity 



In three dimensions we consider the adiabatic collapse of an 
initially non-rotating isothermal gas sphere. This is a stan- 
dard test for SPH codes (Hernquist & Katz 1989; Steinmetz 
& Miiller 1993; Nelson & Papaloizou 1994). To facilitate 
the comparison of our results with those by the above au- 
thors, we adopt the same initial model and the same units 
(M = R = G = 1). The system consists of a 7 = 5/3 gas 
sphere, with an initially isothermal density profile: 



p{r) 



M[R) 1 
2tiR 2 r'' 



(20) 



where M(R) is the total mass inside the sphere of radius R. 
Following Evrard (1988), the profile is obtained stretching 
an initially regular cubic grid by means of the radial trans- 
formation 



r old q 



(21) 



Alternatively it is possible to use the acceptance- 
rejection procedure by Hernquist & Katz (1989). The total 
number of particles used in this simulation is 2176. All the 
particles have the same mass. The initial specific internal 
energy is set to u — 0.05GM/R. 



In PD-SPH, the newtonian forces are calculated using the 
Barnes-Hut tree-code (Barnes & Hut 1986). Tree-code is 
also used to perform the search of neighbors. We always 
include the multipole expansion, and smooth out the accel- 
eration and the potential by means of the Hernquist & Katz 
(1989) formalism. We adopt the opening angle 9 = 0.8, and 
tie the gravitational softening parameter e to the particle 
number by looking at the mean inter-particle separation at 
half- mass radius. In the simulations below e w 0.15. 

4.2 Description of the tests 

The temporal evolution of the system is shown in the various 
panels of Figs. ^| and Q which display the radial velocity and 
the specific internal energy, respectively. Each panel shows 
the variation of the physical quantity under consideration 
(in suitable units) as a function of the normalized radial 
coordinate at different time steps. These (in units of the 
dynamical time scale of the system) are exactly the same as 
in Hernquist & Katz (1989). 

The initial low internal energy is not sufficient to sup- 
port the gas cloud which starts to collapse. Approximately 
after one dynamical time scale a bounce occurs. The system 
can be described as an isothermal core plus an adiabatically 
expanding envelope pushed by the shock wave generated 
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Table 1. CPU time percentage in the adiabatic collapse. 



Section of the code 


Hernquist &: Katz 


This work 


Gravitational computation 


43.28 


42.63 


SPH computation 


27.63 


26.37 


Search of the nearest neighbors 


24.43 


26.22 


Tree construction 


4.50 


4.28 


Miscellaneous 


0.16 


0.50 



at the stage of maximum compression. After about three 
dynamical times the system reaches virial equilibrium with 
total energy equal to a half of the gravitational potential 
energy. 

The temporal evolution of the potential and kinetic en- 
ergies, and the total angular momentum is shown in Fig. || 
As pointed out by Nelson & Papaloizou (1994) the energy 
conservation, expressed as ^ is ^ 8%. This uncertainty 
is caused by neglecting V/i terms and using multiple time 
steps. In brief, when the number of active particle gets too 
small, this causes fluctuations that may affect the energy 
conservation. Evolving the system with an unique time step 
for all the particle lowers the uncertainty on the energy con- 
servation down to about 5 — 6%. As far as the angular mo- 
mentum is concerned, this is conserved within one percent. 
Going back to the panels of Figs. 6 and 7 the present results 
agree fairly well with the mean values of the Hernquist & 
Katz (1989) models provided that the effects of the different 
number of particles are taken into account. For example the 
shock is located at the radial distance 0.18 < r/R < 0.25 
in the models with age t « 0.88 (cf. the velocity panels of 
Fig. ^|). The thermal energy slowly increases for r/R < 0.1 
(cf. the panels of Fig. |^). Finally, at the time of maximum 
compression (t ~ 1.1), the dynamical range in the time steps 
is 20 : 1. The scatter shown by the velocity and internal en- 
ergy in the post shock phases is caused by the relative small 
number of particles, and the use of the multiple time step 
scheme. Similar scatter is present also in Steinmetz & Muller 
(1993) and Nelson & Papaloizou (1994). 
The above numerical tests have been performed on a Digital 
ALPHA OSF/1 workstation. It took about 5400 seconds 
of CPU time for 1019 time steps. It is worth comparing 
the CPU times required by different sections of our code 
with those by Hernquist & Katz (1989) . The comparison is 
summarized in Table 1. It appears that the fraction of CPU 
time spent in key sections of our code agrees with those of 
the Hernquist & Katz (1989) code. 



5 DARK MATTER 

We consider now a collapsing cloud of gas and dark matter 
with initial solid body rotation. The inclusion of dark matter 
is an important step, because it is known to play a key role 
in any modern theory of galaxy formation and evolution (cf. 
Frenk et al. 1996, Persic & Salucci 1997). 

Initially, the dark matter particles have the same spatial 
distribution of the gas particles. Both in fact are obtained 
stretching an initial regular cubic grid (Evrard 1988). There- 
fore, dark matter and gas initially obey the same density 
profile. 

This is a source of a technical difficulty with the tree-code 



algorithm because when two particles occupy the same po- 
sition, the tree-code is not able to distinguish among them, 
when trying to divide the volume into elementary cells each 
of which containing either only one particle or no particles 
at all. This difficulty does not of course occur with the bi- 
nary tree-code algorithm (cf. Benz et al. 1989). To handle 
this problem within a single tree-code scheme, we force the 
tree-code to consider a subdivision as an elementary cell 
also when two particles are enclosed provided that they are 
of different nature. This implies that every particle needs an 
additional flag specifying the species. 

The system under consideration simulates a spherical 
galaxy whit total mass and radius of 1O 12 M0 and 100 kpc, 
respectively. In addition, the system is supposed to be made 
of equal numbers (2176) of particles in form of baryons and 
dark matter but with different individual mass such that the 
total mass fraction of baryons is 0.1 the total mass of the 
galaxy. 

We adopt an initial solid body rotation around the Z- 
axis with angular velocity ui ~ 0.5 which corresponds to a 
frequency of about 1 complete rotation in 10 free-fall time 
scales. So the effects of rotation get sizable during the col- 
lapse. This rotation corresponds to adopting the dimension- 
less spin parameter 

A - GM^ - °-° 8 

is somewhat larger than the value expected from the 
tidal torque in the hierarchically clustering universe theories, 
i.e. A w 0.05 (White 1984, Steinmetz & Bartelmann 1996). 

The gravitational softening parameters e for gas and 
dark matter are in code units (see below) 0.002 and 0.010 to 
which 200 pc and 1 kpc correspond, respectively. They have 
been fixed on the basis of the Evrard (1988) relation 

Grn k< GM 
e R 

The parameter e is kept constant in time for each individual 
particle. However if in a cell two particles of different nature 
(gas and DM) co-exist, the softening parameter is the mass 
weighted mean value. 

Tests similar to those described below can be found in 
Navarro & White (1993) to which our results can be com- 
pared. In this context, the most salient result of the Navarro 
& White (1993) calculations is the different evolutionary be- 
haviour of the dark matter and baryonic components. The 
same results are recovered here as shown in Figs. ^ and JToJ, 
respectively. In our models, however, the difference is less 
pronounced because of the lower initial spin as compared to 
A = 0.10 in Navarro & White (1993). In brief, dark matter 
gets a more flattened spatial distribution than the baryonic 
component. In fact equidensity contours of the baryons must 
coincide with equipotentials which are much rounder than 
contours of equal total density. It is worth recalling that 
observational hints for flatter dark halos have been found 
in galaxies like NGC 5907 (Sackett et al. 1994), NGC 4244 
(Oiling 1996), and NGC 891 (Becquaert & Combes 1997). 



6 NON ADIABATIC PROCESSES 

As far as the evolution of the gaseous component is con- 
cerned, all non adiabatic processes but artificial viscosity 
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Figure 9. Adiabatic collapse: the evolution of the dark matter component as a function of time in the X-Z plane. Time is in code units 



are of paramount importance. They drive in fact the en- 
ergy and momentum loss or gain of these "particles". Since 
the final target is to follow the evolution of a real galaxy, 
processes like cooling and heating, star formation, and feed- 
back play a dominant role, and must be included from the 
very beginning in order to obtain results worth to be com- 
pared with observational data. Unfortunately, some of these 
processes are not yet well understood. The physics of star 
formation, in particular, is far from being assessed so that 
parameterizations of this basic process are imposed. 

In this section we present in some detail our assump- 
tions concerning thermal conduction and cooling processes 
of the gaseous component, which are known to regulate the 
energy budget of a fluid element. They are included in the 
energy equation as source terms. 

6.1 Thermal conduction 

Thermal conduction is calculated according to the formalism 
developed by Monaghan & Lattanzio (1991) 



The translation of it into the SPH language is 



-V ■ (pqVu). 



(23) 



E 1 



(gi - qj)iui - Uj)rjj ■ VW(rij,h) 



(24) 



where q is the thermal conductivity function. This is in turn 
evaluated as 



(25) 



In the above relations, r\ is a parameter securing that the 
denominator of relation (25) remains different from zero, 
and and g is another suitable parameter set equal to 0.25. 

Thermal conduction is expected to be effective in sit- 
uations of strong shocks, and/or "wall heating" when two 
flows or two fluid elements collide. Thermal conduction is a 
source of heating in the energy equation (see below) . 



6.2 Cooling 

Radiative cooling is the crucial mechanism responsible for 
condensation and collapse of baryonic gas into galaxies and, 
inside galaxies, for the occurrence of star formation. Cooling 
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Figure 10. Adiabatic collapse: the evolution of the baryonic component as a function time in the X-Z plane. Time is in code units 



processes are rather well known, and have already been in- 
cluded in current models of galaxy formation and evolution 
albeit in different detail. See for instance Katz, Weinberg 
& Hernquist (1995) for a very accurate treatment of cool- 
ing processes of different nature. Among the various cooling 
agents we recall: (i) the radiative cooling by atomic and 
molecular processes; (ii) the inverse Compton mechanism in 
presence of microwave background radiation. 

Radiative cooling. The most commonly used radia- 
tive cooling functions are those by Katz & Gunn (1991) 
because they are analytical and continuous with their first 
derivatives. These radiative cooling functions are calculated 
assuming the typical helium abundance by mass Y=0.25. 
They do not contain, however, the dependence on metal- 
licity which in contrast is expected to be important dur- 
ing the evolution of a galaxy, as a natural consequence of 
the chemical enrichment. To cope with this drawback of 
the Katz & Gunn (1991) cooling functions, we adopt here 
those elaborated by Chiosi et al. (1997). In brief, these cool- 
ing function Ac(«i, pi, Zi) are derived from literature data, 
and include different radiative processes. For temperatures 



greater than 10 4 K they lean on the Sutherland & Dopita 
(1993) tabulations for a plasma under equilibrium condi- 
tions and metal abundances [Fe/H]=-10 (no metals), -3, -2, 
-1.5, -1, -0.5, (solar), and 0.5. For temperatures in the 
range 100 < T < 10 4 the dominant source of cooling is 
the H2 molecule getting rotationally or vibrationally excited 
through a collision with an H atom or another H2 molecule 
and decaying through radiative emission. The data in use 
have been derived from the analytical expressions of Hollen- 
bach & McKee (1979) and Tegmark et al. (1996). Finally, 
for temperatures lower than 100 K, starting from the re- 
lation of Theis et al. (1992) and Caimmi & Secco (1986), 
they corporate the results of Hollenbach & McKee (1979), 
and Hollenbach (1988) for CO as the dominant coolant. The 
following analytical relation in which the mean fractionary 
abundance of CO is given as a function of [Fe/H], is found to 
fairly represent the normalized cooling rate (i.e. Ac (CO) /n 2 
with n the number density of particles) 

AciCO) = i.exio-^io^^l- 1 - 699 ^ - 5 erg cm 3 B -\(26) 

The normalized cooling rate (in units of erg cm 3 s _1 ) over 
the whole temperature range is shown in Fig. O, in which 
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Figure 11. Collapse of a clump of gas and dark matter. Upper panel: conservation of the potential (dashcd-dottcd line), kinetic (dotted 
line), thermal (dashed line), and total energy (continuous line). Lower panel: conservation of the total angular momentum 



the effects of the metallicity are clearly visible. It is worth 
mentioning that no re-scaling of the cooling rate from the 
various sources has been applied to get the smooth curves 
shown in Fig. |l2| For the purposes of comparison, we display 
in Fig. 11 the cooling rate of Katz & Gunn (1991). The ma- 
jor point of disagreement is soon evident because their cool- 
ing rate which by definition is for a zero-metal composition 
actually lies well above the zero-metal curve of Sutherland 
& Dopita (1993). This implies that Katz' & Gunn's (1991) 
cooling is unphysically more efficient in situations of nearly 
zero metallicity, such as the protogalactic phase. 

Inverse Compton. We do not take into account in- 
verse Compton cooling because of its strong dependence on 
redshift (Ikeuchi & Ostriker 1986). During the protogalac- 
tic phase of the evolution of a galaxy, from which we start 
our simulations, the inverse Compton cooling has already 
got much lower than the radiative cooling, so that it can be 
neglected. 

Although the cooling rate has been derived over a wide 
range of temperatures from a few degree to above 10 8 , in 
reality only the portion above a suitable temperature T min 
can be used in model calculations because of the maximum 
resolution achievable in the N-body simulations. This limit 
is derived from imposing that the Jeans mass does not fall 
below a critical value, which is conventionally assumed to be 



the mass of four gas particles (cf. Katz & Gunn 1991). The 
limit temperature is given by the Bonnor-Ebert expression 
(Ebert 1955, Bonnor 1956) 



= 0.89553 



4m ? 
Pa 



2/3 



2/3' 



P a , (27) 



where m p is the mass of a baryonic particle, k is the Boltz- 
mann constant, [i and e are the mean molecular weight and 
the gas gravitational softening parameter, respectively. 

The limit temperature T m i n is locally computed for ev- 
ery gas particle. Finally, the cooling functions are stored in 
the code as look-up tables. 



6.3 Heating 

Heating of the gas is caused by many processes among which 
we consider the thermalisation of the energy deposit by su- 
pernova explosions (both Type I and II), stellar winds from 
massive stars, the ultraviolet flux from massive stars, the 
cosmic ultraviolet flux, and finally sources of mechanical na- 
ture. 

Supernovae. The rate of supernova explosions 
RsNi,n(t) over the time interval At is calculated according 
to standard prescriptions (cf. Bressan et al. 1994; Tantalo 
et al. 1996, 1997; Greggio & Renzini 1983 for Type I SN 
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Figure 12. Net cooling rate Ac/n 2 in ergs cm 3 s 1 as a function of temperature, n is the number density of particles. The cooling rate 
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& Secco (1986) and Theis et al. (1992) for T < 100 K. Each curve corresponds to a different metallicity as indicated. For comparison the 
cooling rate by Katz & Gunn (1991) is superposed (open triangles) 



in particular). Knowing the amount of energy released by 
each SN explosion, the total energy injection over the time 
interval At is the sum of two terms, Esni and Esnii, of 
type 

Esni, ii = / esNi,iiRsNi,n (t')dt', (28) 

where RsNi,n{t) is the number of supernovae per unit time, 
and esNiji is the energy injected per SN explosion. This 
formulation slightly differ from the one commonly adopted 
in chemical models of galaxies, in which esNi,n(t) is the 
thermalisation law of SN remnants include of the cooling 
processes (cf. Bressan et al. 1994, Tantalo et al. 1996, and 
Gibson 1995). In our scheme, first cooling effects are left to 
the energy budget equation, second owing the limitations 
implicit in the N-body technique we are not yet able to re- 
solve the star forming fluid to the level of individual stars but 
only to that of massive cluster-like structures in which star 
can be formed according to a given initial mass function (see 
below). Therefore, following in detail the thermalisation law 
of SN is impossible. A reasonable compromise is obtained 
by using sufficiently small time steps (cf. also Chiosi et al. 
1997 for a similar approach). The explicit formulations for 
RsNi,n(t) will be presented in the next section. 

Stellar winds. The rate of energy injection by stellar 
winds is 

E w = I ewRw{t')dt', (29) 

J At 



where Rw is the number of stars per unit time expelling 
their envelopes during the time interval At and, in analogy 
with the SN remnants, ew is the kinetic energy of stellar 
winds. A losing mass star is expected to deposit into the 
interstellar medium the energy 

*o = ?x»irx,(M) J 1 (30) 

2 Zq 

where M e j(M) is the amount of mass ejected by each star 
of mass M, v(M) is the velocity of the ejected material, and 
r\ is an efficiency factor of the order of 0.3 (Gibson 1994 
and references therein). The calculation of the rate Rw(t) 
is postponed to the next section. 

Ultraviolet flux from massive stars. The rate of 
energy injection from the ultraviolet flux emitted by massive 
stars is 

E vv = / e uv Ruv(t')dt', (31) 

J At 

where Ruv is the number of massive stars per unit time 
whose mass is the range 10 -f- 120 Mq, and euv is the 
amount of ultraviolet energy emitted by each star. To cal- 
culate euv the following procedure is adopted. We suppose 
that massive stars are located on the zero age main sequence, 
i.e. obeying well known mass-luminosity-effective temper- 
ature relationships, L/L Q (M/M Q ) and T e f f (M/M @ ). In 
principle, there should be an additional dependence on 
the chemical composition which, however, can be neglected 
here. The relationships L/L Q (M/M Q ) and T eff (M/M Q ) 
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Figure 13. Testing cooling: the final distribution in X — Z plane for gas particles in simulations with different cooling prescriptions. 
The right panel is for the cooling rate by Katz & Gunn (1991), whereas the left panel is for the cooling prescription described in the text 



are taken from the library of stellar models/isochrones of 
Bertelli et al. (1994 and references therein). For the sake of 
simplicity we can approximate the spectral energy distribu- 
tion of any such stars as pure black body emission of given 
T e ff, for which we can immediately estimate the fraction 
Fuv(M) of flux emitted short- ward of 4000 A (our range 
for UV light) by a star of mass M. The UV flux emitted by 
each star is 



Wv = Fuv (M) x L(M), 



(32) 



where L is the total luminosity of the star. Once again the 
calculation of Ruv{t) is postponed to the next section. 

Cosmic UV radiation. As far as UV radiation field 
is concerned, its effects have been thoroughly investigated 
by Navarro & Steinmetz (1996) who reached the conclusion 
that the cosmic UV radiation plays some role only when 
the gas has already accreted onto the protogalaxy from the 
surrounding medium. Heating by cosmic UV radiation is not 
included in the present models. In this context, it is worth 
recalling that all models in which the cosmic UV radiation 



is the sole source of heating face the so-called overcooling 
problem (cf. Navarro & Steinmetz 1996). 

Total radiative heating. The total heating rate due 
to radiative processes Hr(ih, pi, Zi) to be included in the 
energy budget equation (see below) is 



H R = 



Esni + Esnii + Ew + Euv 
At 



(33) 



Mechanical heating. In principle there are various 
mechanical processes that contribute to heat gas particles. 
Two of them have been explicitly taken into account in the 
basic energy equation (^). They have been shortly referred 
to as Ti t M and are general functions of type Vi t M(ui, pi, Zi). 



6.4 More details on the energy equation 

At this stage of the analysis we can write down the detailed 
expression of the energy equation in all its components. To 
this aim we need only to convert the radiative heating into 
the conventional SPH formalism 
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Figure 14. Cooling radiation (in code units): the solid line is with the cooling rate by Katz & Gunn (1991), whereas the dashed line is 
for the cooling rate used in this paper 



Twiui, pi, Zi) =y^ l H Rd x ^W(i,j) (34) 
i 

with obvious meaning of the symbols. The final energy equa- 
tions is 

dui A c ,i(ui,Pi,Zi) . . 

—TT = 1 M,i + 1 R,i(Ui,pi, Zi) (6b) 

at pi 

with 

^R,i(ui,pi,Zi) — Hc,i( u i, Pu Zi) + H R ,i(ui, pi, Zi) (36) 

where Yr^ the sum of thermal conductivity and all radiative 
heating sources. The energy m is per unit mass, whereas all 
other quantities are per unit mass and time. 

The functional dependence of the cooling term imposes 
an implicit solution of the energy equation. This is per- 
formed using an hybrid scheme, i.e. a combination of the 
Newton-Raphson and bisection methods. The associated 
second order updating of the thermal energy follows the 
technique of Hernquist & Katz (1989). 

Finally, it is worth mentioning that under high cooling 
efficiency situations can be met in which the energy may be- 



come negative. To cope with this unphysical result of mere 
numerical nature we impose that a gas particle cannot lose 
more than a half of its thermal energy per time step (Hern- 
quist & Katz 1989). 



7 TESTING COOLING PRESCRIPTIONS 

In this section we show an ideal galactic model in which 
no star formation, no chemical enrichment, and no energy 
deposit from the various sources are let occur. This model 
is meant to show the net results of different prescriptions 
for the cooling rate. The model has zero-metal content and 
apart from the cooling rate, it is identical to that by Katz 
& Gunn (1991) so that the comparison is possible. 

We simulate the collapse of a spherical galaxy with 
mass of 10 12 Mq and radius of 100 kpc. Like in the pre- 
vious test the baryon fraction is 0.1, while the initial gas 
temperature is 10 4 °K. In this simulation, time, spatial co- 
ordinates, energy, and density are expressed in the following 
units: [time] = 7.41 x 10 8 year, [space] = 100 kpc, [energy] = 
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Figure 15. Conservation of energy and angular momentum in a collapsing galaxy with the cooling rate adopted in this paper. The 
model is the same as in the left panel of Fig. 13 



"erg gr _1 , and [density] = 1.61 x 10 _25 gr cm -3 . 
The temporal evolution of the models has been followed up 
to the age of 2.2 dynamical times ( approxl.&Gyr). 

The results are schematically presented in Figs, [l^ and 
K3 . Fig. [l3] shows the distribution of particles on the X — Z 
plane at three different ages (in code units) for our cooling 
rate, left panels, and Katz & Gunn (1991), right panels, 
whereas Fig. [l4| compares the cooling radiation for the two 
assumptions for the cooling rate. 

As the zero-metal cooling rate of Sutherland & Dopita 
(1993), on which our prescription stands for temperatures 
above 10 4 K, is much less efficient than that of Katz & 
Gunn (1991), this easily explains why significantly longer 
time scales are need to form disk structures (cf. Fig. 
and the much lower amount of radiated energy (cf. Fig. 

Finally, in Fig. jHs] we show the conservation of total 
energy and angular momentum as a function of time for 
the case with the new cooling rate. The small decrease of 
the total energy at increasing time is caused by the lost 
of thermal energy radiated away by cooling processes (for 
comparison see the detailed trend of energies in Fig. 11). 



Angular momentum is conserved within 5% accuracy (see 
also Fig. 11 for comparison). 



8 GOING TOWARD COMPLETE MODELS: 
STAR FORMATION, FEED-BACK, AND 
CHEMICAL ENRICHMENT 

Implementing star formation into N-body simulations of 
galaxies is a cumbersome affair, reflecting partly our poor 
knowledge of this fundamental process, and partly the lim- 
itations still imposed by N-body technique itself. 

To accomplish the task two basic steps are required: 
firstly we need to determine the star formation rate (SFR), 
and secondly we must treat the effects of the star formation 
on the interstellar medium (the so-called feed-back). 



8.1 Star formation rate 

To include star formation in the N-body scheme two impor- 
tant informations are needed: firstly we have to set suitable 
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Figure 16. Simulation of a disk-like galaxy with no feed-back: the spatial distribution of gas (left panel) dark matter (central panel) 
and stars (right panel). Starting from the bottom the three panels of each raw refer to ages of 0.5, 2.5, and 5 Gyr 



conditions under which stars can form, secondly we must 
express the SFR law governing the efficiency at which gas 
is turned into stars (cf. Steinmetz 1995 for a review on the 
subject). 

When to form stars? There are at least three physical 
conditions to be met in order to activate the star forming 
process. First, a fluid element is eligible to make stars if it is 
part of a convergent flow, i.e. the particle velocity divergence 
must verify the condition 

V • v, < 0. (37) 
Second, the fluid element must be Jeans unstable 

hi 



pi per iti 



(40) 



> 



Finally, the gas particle must remain cool 

Tcool < Tff. 



(38) 



(39) 



This last condition is usually translated into an overdensity 
criterion. For instance, Navarro & White (1993) uses the 
following expression 



where p cr it — 7 X 10~ g cm . This minimum threshold 
density actually holds only in the case of a single cooling 
function (no dependence on chemical composition). How- 
ever, as star formation causes metal enrichment of the in- 
terstellar medium, the threshold moves towards lower values 
at increasing metallicity. Because of this, we prefer to adopt 
condition (B9) which is of general validity. 

Star formation law. A very popular, empirically 
based, rate of star formation is the Schmidt (1959) law, ac- 
cording to which 



*(*) 



dpgas 

dt 



dp s 



dt 



Pgas 
T ff 



(41) 



where c* is the specific efficiency (in most cases a free pa- 
rameter) . 

In dense regions (the ones prone to star formation), the 

cooling time t coo i is typically much shorter than the dynam- 

r /2 

ical time Tff. Since Tff oc p gaa , the SFR is proportional to 
p^as, a result that agrees with Schmidt's (1959) observa- 
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Figure 17. Simulation of a disk-like galaxy with no feed-back: the same as in Fig. 16 but onto the X-Y plane 



tional estimates. In most numerical simulations, the specific 
efficiency is taken to be c* = 0.1. 



8.2 Star-like particles 

When star formation has started, at any time-step At a 
star-like particle is created, whose mass is 



x (1 — exp(- 



c+At 

Tff 



))■ 



(42) 



and the mass of the parent gas particle is consequently re- 
duced. The above relation simply follows from integrating 
equation (|ll|) over the time step At. The mass of the star- 
like particles obviously depends on the resolution (number 
of particles) of the simulation. 

The star-like particle is supposed not to immediately 
acquire its own individuality, but to leave the parent gas 
particle only when the mass of this latter falls below a cer- 
tain value (typically 50% of the original mass) . Furthermore, 
recurrent episodes of star formation within the same gas par- 
ticle are possible, so that gas is depleted in discrete steps. 
Finally, the star-like particle it treated as a collisionless ob- 



ject. Therefore, the element of fluid in which star formation 
occurs is considered as a hybrid particle, whose collisional 
component experiences both hydrodynamical and gravita- 
tional forces, while its collisionless component feels only the 
gravitational field. 

To avoid non physical accelerations due to the sudden 
decrease of the gas mass, particles in which star formation 
is active are stored in the lowest time-bin. 

When a star-like particle leaves its parent, it inherits 
its position in the phase space, and its gravitational soft- 
ening parameter. Resolution limitations force on us to con- 
sider only four star formation events per gas particle, like in 
Navarro & White (1993). 



8.3 Feed-back 

To quantitatively evaluate the feed-back by all heating pro- 
cesses already outlined in the previous section, we need to 
know the initial mass function, according to which the newly 
borne stars will distribute by number in each mass interval 
dM. For the sake of simplicity we adopt here the classical 
Salpeter (1955) law 
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Figure 18. Simulation of a disk-like galaxy with no feed-back: 
the rates of star formation (top panel) and supernova explosions 
(Type II middle panel, Type I bottom panel) 



$(M)dM = A ■ M~ x dM, 



(43) 



where x = 2.35 (the Salpeter value) and A is the normal- 
ization constant. This is derived from imposing that stars 
can be formed with masses in the range 0.1 4- 12OM0, and 
assuming that the integral of eq. ( fl3| ) over this mass range 
is equal to one. It is worth recalling that only stars more 
massive than about 0.9 Mq will be able to pollute the in- 
terstellar medium over a time scale shorter than the Hubble 
time (cf. Tinsley 1980, Matteucci 1997). The normalization 
constant is A = 0.06. Other IMFs are already implemented 
into the code (Miller & Scalo 1979). 

Supernova rates In the standard scenario of stellar 
evolution (cf. Woosley 1986, Chiosi 1986; Greggio and & 
Renzini 1983) Type II supernovae occur in single stars more 
massive than say 8 Mq, whereas Type I (more precisely la) 
supernovae occur in binary stars containing a re-juvenated 
white dwarf brought to explosion via the C-deflagration 



1= 

CD 




Figure 19. Simulation of a disk-like galaxy with no feed-back: 
the age-metallicity relation for the gas component in the galactic 
disk. Filled circles show the age-metallicity relation of Disk stars 
in the Solar Vicinity by Edvardsson et al. (1993) 



mechanism by mass accretion from the companion. Further- 
more, Type II supernovae from progenitors more massive 
than say 30 Mq leave behind a black hole (scarcely con- 
tributing to the enrichment of the interstellar medium in 
heavy elements), whereas Type II supernovae from progen- 
itors in mass range 8 to 30 Mq leave a neutron star and 
most effectively contribute to chemical enrichment in heavy 
elements. Type I supernovae have binary system progenitors 
with typical total mass in the range 3 to 16 Mq. 

To evaluate the rate of Type I occurrence one has to 
know the percentage of binary systems with respect to single 
stars. Let C be this fraction. Following Greggio & Renzini 
(1983), the rate of Type I supernovae is 



1/2 



C 



■yr+l 



$ B (M B )dM B x 



(l + 7 ) M 7 *[i-t (flMB) ]d it . 



(44) 



In the above expression Mb is the total mass of the binary 
system, $b(Mb) is the IMF of binary stars, which is as- 
sumed here to be the same as for single stars [eq.(^)] for 
the sake of simplicity, jx = is the fractionary mass of 
the secondary, which drives the evolution, and 7 is a coef- 
ficient, usually equal to 2. Finally, M B>in f and Mb, max are 
the lower and upper limits for the total mass of the binary 
stars, respectively, and fj,i„f is the lower limit for the frac- 
tionary mass of the secondary. These masses are calculated 
with the following receipt 



M B jn f = MAX(2M 2 ,2), 

M B ,max =8 + M 2 , 



(45) 
(46) 
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Figure 20. Simulation of a disk-like galaxy with full feed-back: the spatial distribution of gas (left panel) dark matter (central panel) 
and stars (right panel). Starting from the bottom, the three panels of each raw refer to ages of 0.5, 2.5, and 5 Gyr 



IH nJ = MAX(M 2 (t)/M B , 1/2). 



(47) 



Finally, <3/(i) is the SFR at the time equal to the difference 
between the current time and the lifetime of Mi. 



The total rate of Type II supernovae is 

(•120 

Rsnii = / <S>{M)Psi(t M )dM+ 

/•M Bm „ 

(1-C) / $(M)*(i M )dM 



(48) 



with obvious meaning of the symbols. 



Finally, the constant C above is calibrated looking at 
the rate of Type la supernovae in spiral galaxies (van den 
Bergh & McClure 1994). C turns out to be w 0.04048. The 
contribution by supernovae explosions to feed-back has been 
also included in the Tree-SPH model of the Milky Way by 
Raiteri et al. (1996). 

Stellar winds and UV fluxes. The rates Rw and 
Ruv, i.e. number of stars per unit mass and time contribut- 



ing to the energy injection by stellar winds and UV flux are 
simply given by 



Ri 



<&{M)^(t M )dM. 



(49) 



All mass-lifetime relationships entering the above cal- 
culations of the rates RsNiji, Rw, and Ruv are derived 
from the library of stellar models/isochrones of Bertelli et 
al. (1994). Whenever possible the effect of the chemical com- 
position of the fluid element is taken into account. 

Final remark. Each supernova explosion produces 
10 51 ergs of energy which is injected in the interstellar 
medium in form of kinetic energy. A comparable amount 
of energy is supplied by a massive star during its whole life- 
time in form of stellar wind. In the case of supernovae, only 
a small fraction of this energy is thermal (a few percent), 
and this small fraction is immediately (100 years is the typ- 
ical time scale) radiated away for instance in form of X-rays 
emission. The effect of this energy input on the velocity field 
of the surrounding medium is expected to negligible. In brief, 
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Figure 21. Simulation of a disk-like galaxy with full feed-back: the same as in Fig. 20 but onto the X-Y plane 



while the typical space resolution of current simulations is 
about 1 kpc or more, supernova shocks may affect volumes 
of about 100 pc radius around the exploding star (Mazzali 
& Cappellaro, private communication) Therefore it is not 
possible to follow the dynamics of such explosions down to 
the required resolution. On the base of these considerations, 
the claim that some fraction of the supernova energy budget 
(10~ 4 or so of the total) can modify the velocity field of the 
neighboring volumes of our star-like particles is not solidly 
grounded. We prefer to simply assume that the cumulative 
effect of the supernova explosions increases the mean tem- 
perature and internal energy of the gas on a scale larger than 
the size a fluid element. 



8.4 Chemical evolution 

Following , in our simulations, the chemical enrichment of 
the inter-stellar medium follows the prescriptions by Stein- 
metz & Miiller (1994). It stands on the closed-box, instanta- 
neous recycling description (cf. Tinsley 1980), i.e. the metal 
abundance Z of a fluid element is computed as 



AZ t = -Y z Am ^ gas . 



(50) 



where Yz is the so-called yield per stellar generation. Given a 
certain stellar nucleosynthesis scenario (cf. Matteucci 1997), 
Yz can be calculated a priori in a way mutually consistent 
with the assumed IMF. Using the results of chemical models 
by Portinari et al. (1997) a good assumption for yield per 
stellar generation is Yz = 0.005. 

Once metals are synthesized, they are spread according 
to the SPH formalism: 



<z i >=y Zj ^. 

r-1 Pj 



W(Ra,h). 



(51) 



Groom(1997) suggests that the chemical enrichment of 
the interstellar medium is better described by the diffusive 
equation 

dZ 



dt 



-kV Z 



(52) 



where Z is the metallicity, and n the diffusion coefficient, 
with k equal to 50 km sec -1 60 pc. The SPH translation, 
which requires second derivatives of the kernel, is 
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aZi sr^ 



3 = 1 



Pj + Pi 



IV WijQfi - rjl.hi. 



(53) 



Chemical enrichment is expected to occur on the same 
spatial scale of the energy release of supernova explosions, 
i.e. about 100 pc, which is much smaller than the spatial 
resolution of our model. 

Although the use of 2 nd derivative of the smoothing 
function may introduce spurious fluctuation in the resulting 
metallicity implicit in equation (55), this is less of a problem 
as far as the global enrichment in metals is concerned. 

In fact, in semi-analytical models of chemical evolu- 
tion of disk-like galaxies (Tinsley 1980, Carraro et al. 1997, 
Portinari et al 1997, and references therein) the metallicity 
is know to quickly reach the metallicity yield of the con- 
tributing stellar population. 

Therefore we expect the result not to critically depend 
on the particular scheme used to evaluate the metallicity 
variation with the SPH method. This is implicitly confirmed 
by the agreement between theoretical results and observa- 
tional data for the solar vicinity that are presented in the 
section below. 

We have adopted equation (55), and converted the 
metallicity Z to [Fe/H] by means of the relation 

[Fe/H] = log(Z) + 1.739 (54) 

taken from Bertelli et al. (1994). 

9 THE CASE OF A DISK-LIKE GALAXY 

The code described in the previous section has been used to 
study the formation and evolution of disk-like galaxies. The 
initial model consists of a 10 12 Mq spherical galaxy enclosed 
within a radius R = 100 kpc and containing 1000 particles of 
baryonic mass and 1000 particles of dark matter. We have 
assumed that the fractionary total mass of of baryons is 
0.1 the total mass of the galaxy and the spin parameter 
A = 0.08, which is in the range of current estimates for disk 
galaxies. 
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Figure 22. Simulation of a disk-like galaxy with full feed-back: 
the rates of star formation (top panel) and supernova explosions 
(Type II middle panel, Type I bottom panel) 



9.1 No feed-back 

In Figs. |l6| and ^ we show the dynamical evolution of a 
disk-like galaxy in absence of feed-back. Under the combined 
action of rotation and cooling, gas settles onto a disk (cf. 
left panels of Figs. |l6| and [l?]), and stars start to be formed. 
Remarkably, at the age of about 5 Gyr a bulge-like structure 
made of stars is observed, whereas a disk-like structure is 
not yet visible because of the insufficient resolution, and our 
assumption that stars are considered as such when the gas 
fraction of each star forming unit has fallen below 0.5 the 
total mass. We expect the disk to appear at later ages. Like 
the case of the collapsing galaxy without cooling and star 
formation, a flattened halo of dark matter can develop (cf. 
the central panels of Figs [li] and [It] and Fig.[]). Another 
important result of the new cooling rate is that gas particles 
survive in the halo for significantly long periods of time. 

The temporal evolution of the star formation and su- 
pernova rates is shown in Fig. As expected the rate of 



star formation is low and nearly constant with time. Fur- 
thermore, while Type II supernovae closely follow the trend 
of the star formation, Type I SN appear in significant pro- 
portions only later than about 1.5 Gyr due to the rather long 
lifetime of their progenitors. Finally, in Figjl^, we show the 
age-metallicity relation for the gas component on the disk 
of the galaxy and compare it with the age-metallicity re- 
lation for disk stars in the solar vicinity by Edvardsson et 
al. (1993). Despite the crudeness of our modeling of chemi- 
cal enrichment the agreement is remarkable. In particular at 
ages older than 2 Gyr, when the metallicity has saturated 
to the yield. 

9.2 Complete feed-back 

The spatial evolution of the disk-like galaxy in presence of 
full feed-back in shown in Figs. ^ and|5l], which display the 
gas (left panels), dark matter (central panels), and star par- 
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Figure 23. Simulation of a disk-like galaxy with full feed-back: Figure 24. The gas fraction in models with no feed-back (dashed- 

the age-metallicity relation for the gas component on the galactic line) and models with full fee-back (solid line) 
disk. Filled circles show the age-metallicity relation for Disk stars 
in the Solar Vicinity by Edvardsson et al. (1993) 



tides (right panels) at three different ages, namely 0.5 (bot- 
tom), 2.5 (middle) and 5 Gyr (top). As expected the galaxy 
remains less spatially concentrated owing to the energy in- 
put from the various sources and consequent heating. This 
is particularly true for the gas and star particles, whereas 
no sizable effect is seen on the dark matter component. The 
rates of star formation and Type I and II supernova explo- 
sions as a function of time are shown in the three panels of 
Fig. Although the global trend is as in the previous case, 
now the three rates are down by approximately a factor of 
two. This can be easily understood as a result of heating 
which keeps the system less spatially concentrated (lower 
density) thus immediately lowering the rate of star forma- 
tion and supernova explosions in turn. 

In Fig.E3l we show the age-metallicity relation for the 
gas component on the disk of the galaxy and compare it 
with the age-metallicity relation for disk stars in the solar 
vicinity by Edvardsson et al. (1993). The same remark made 
on Fig. 19 still applies. 

Finally, in order to prove that the mean density of the 
two model galaxies is actually different we plot in Fig. |24| 
the ratio m gas /(m gas + m sta rs) as a function of the age. The 
dashed line is the galaxy with no feed-back at all, whereas 
the solid line is for the galaxy with full feed-back. 



10 CONCLUSIONS AND FUTURE 
PERSPECTIVES 

In this paper, we have presented a new Tree-SPH code suited 
to studies on the formation and evolution of galaxies. 

Particular attention has been paid to include modern 
descriptions of non adiabatic processes, such as cooling and 



heating by supernova explosions, and stellar winds and UV 
radiation from massive stars, so that feed-back in the en- 
ergy balance is not a parameter but naturally follows from 
the assumed star formation rate and initial mass function. 
Furthermore we have followed the chemical enrichment of 
gas as the consequence of star formation and include the 
dependence of the cooling rate on the gas chemical compo- 
sition (metallicity). These are points of major difference and 
improvement with respect to similar codes in literature. 

Although the initial conditions are not tighten to any 
specific cosmological scenario, still the ones we have adopted 
lead to reasonable results. This point will be the subject of 
future improvements. 

In addition to many classical tests aimed at checking 
the performance of the code and its response to different 
physical assumptions, we have presented the results for two 
disk-like galaxies with and without feed-back for purposes 
of comparison both internal and with similar studies in lit- 
erature. 

Particularly relevant is the role played by the depen- 
dence of the cooling rate on the chemical composition 
(metallicity), which owing to its strong impact on the final 
results can no longer be neglected in simulations of galaxy 
formation and evolution. 

Work is in progress to extend these models to the case of 
elliptical galaxies and to build the proper interface between 
this code and the spectro-photometric code of Bressan et 
al. (1994) so that the spatial history of integrated spectra, 
magnitudes, colors, line strength indices etc. can be derived 
together with the dynamical history for galaxies of different 
morphological type. 
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